##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/DriverChoose/GeneLOH"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_info <- opt$sample_info
out_path <- opt$out_path
ccf_file <- opt$ccf_file

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_info <- data.frame(fread(sample_info))
dat_ccf <- data.frame(fread(ccf_file))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
gene_list <- c("TP53" , "APC" , "PIK3CA" , "CDH1")

###########################################################################################
## 加载分类
dat_ccf <- merge( dat_ccf , unique(dat_info[,c("ID" , "Type")]) , by = "ID" )
dat_ccf2 <- dat_ccf
dat_ccf2$Type <- "All"

dat_ccf <- rbind( dat_ccf , dat_ccf2 )

###########################################################################################
## 判断基因多少比例出现在Trunk、Pre_Private、Inv_Private
## 计算p值
trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

result_plot <- c()

for(typeN in unique(dat_ccf$Type)){
    print(typeN)
    for( geneN in gene_list ){
        print(geneN)
        tmp <- subset( dat_ccf , Variant_Classification %in% Variant_Type & Hugo_Symbol == geneN & Class != "IM" & Type == typeN )

        if(nrow(tmp) > 0 ){
            ## 判断拷贝数改变
            tmp$CNV_type <- "Nochange"
            tmp$CNV_type <- ifelse( tmp$minor_cn == 0 & tmp$total_cn==2 , "LOH" , tmp$CNV_type )
            tmp$CNV_type <- ifelse( tmp$total_cn > 2 , "Gain" , tmp$CNV_type )
            tmp$CNV_type <- ifelse( tmp$total_cn == 1 , "Loss" , tmp$CNV_type )

            res_tmp <- c()
            for( sample in unique(tmp$ID) ) {
                tmp_use <- subset( tmp , ID == sample )
                tmp_use <- unique(tmp_use[,c("ID" , "Share" , "CNV_type")])

                ## 一个人若发生多个突变，算共享的的
                if( length(which(tmp_use$Share==TRUE)) > 0 ){
                    tmp_use <- subset( tmp_use , Share == TRUE )
                }

                ## 三个样本存在多个患者拷贝数改变不一致，选择有改变的样本
                if( nrow(tmp_use) > 1 ){
                    tmp_use <- tmp_use[2,]
                }

                res_tmp <- rbind(res_tmp , tmp_use)
            }

            res_tmp$CNV_change <- ifelse( res_tmp$CNV_type == "Nochange" , "Nochange" , "DoubleHit" )
            res_tmp2 <- data.frame(table(res_tmp$Share , res_tmp$CNV_type , res_tmp$CNV_change))
            colnames(res_tmp2)[1:4] <- c("Trunk" , "CNV_type" , "CNV_change" , "count")
            
            tmp1 <- table(res_tmp$Share , res_tmp$CNV_type)
            if(ncol(tmp1) > 1 & nrow(tmp1) > 1){
                p <- fisher.test(tmp1)$p.value
            }else{
                p <- 1
            }

            if( p < 0.001 ){
                p_text <- trans(p)
            }else{
                p_text <- paste0( "P == " , round(as.numeric(p) , 2) ) 
            }

            res_tmp2$p <- ""
            res_tmp2$p[1] <- p
            res_tmp2$p_text <- ""
            res_tmp2$p_text[1] <- p_text

            ## 计算每一类的样本数量
            res_tmp2 <- res_tmp2 %>%
            group_by(Trunk) %>%
            mutate(count_all=sum(count),
               ratio=count/count_all)

            res_tmp2$Gene_Symbol <- geneN
            res_tmp2$Type <- typeN

            result_plot <- rbind( result_plot , res_tmp2 )
        }
    }
}

out_name <- paste0(out_path , "/DoubleHit.maintained.tsv")
write.table( result_plot , out_name , row.names = F , sep = "\t" , quote = F )

###########################################################################################
## 基因分布堆叠图
result <- data.frame(result_plot)
result$ratio[is.na(result$ratio)] <- 0
result$value_text <- round(result$ratio , 2) * 100

result$CNV_type <- factor( result$CNV_type , levels = c("Gain" , "Loss" , "LOH" , "Nochange") )
result$TrunkUse <- ifelse( result$Trunk == "TRUE" , "Trunk" , "Private" )
result$Gene_Symbol <- factor(result$Gene_Symbol , levels = gene_list )

col <- c(
    rgb(190,97,99,alpha=255,maxColorValue=255),
    rgb(91,139,101,alpha=255,maxColorValue=255),
    rgb(76,123,161,alpha=255,maxColorValue=255),
    "grey"
  )

names(col) <- c("Gain" , "LOH" , "Loss" , "Nochange")
result$value_text <- ifelse( result$ratio == 0 , "" , result$value_text )
result$TrunkUse <- paste0(result$TrunkUse , "(" , result$count_all , ")")
result$TrunkUse <- factor( result$TrunkUse , levels = unique(result$TrunkUse )[order(unique(result$TrunkUse ) , decreasing=T )] )

for(typeN in unique(result$Type)){
    plot <- ggplot( data = subset(result , Type == typeN) , aes( x = TrunkUse , y = ratio , fill = CNV_type ))+
        geom_bar(position = "stack", stat = "identity") + 
        theme_bw()+
        labs(x="",y="Proportion (%)")+
        facet_grid(.~Gene_Symbol,space='free_x',scales='free_x') +
        theme(panel.grid = element_blank())+
        scale_fill_manual(values=col) +
        ylim(0,1.05)+
        geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
        geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                    legend.position ='right',
                    legend.title = element_blank() ,
                    panel.grid.major=element_line(colour=NA),
                    plot.title = element_text(size = 12,color="black",face='bold'),
                    legend.text = element_text(size = 8,color="black",face='bold'),
                    axis.text.y = element_text(size = 7,color="black",face='bold'),
                    axis.title.x = element_text(size = 12,color="black",face='bold'),
                    axis.title.y = element_text(size = 12,color="black",face='bold'),
                    strip.text.x = element_text(size = 7 , face = 'bold'),
                    axis.text.x = element_text(size = 8,color="black",face='bold') ,
                    axis.line = element_line(size = 0.5))

    out_name <- paste0(out_path , "/DoubleHit.maintained.",typeN,".pdf")
    ggsave( out_name , plot , width = 7 , height = 4 )
}

## 提取TP53看拷贝数改变
plot <- ggplot( data = subset(result , Gene_Symbol == "TP53") , aes( x = TrunkUse , y = ratio , fill = CNV_type ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion (%)")+
    facet_grid(.~Type,space='free_x',scales='free_x') +
    theme(panel.grid = element_blank())+
    scale_fill_manual(values=col) +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 8,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold') ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/DoubleHit.maintained.TP53.pdf")
ggsave( out_name , plot , width = 6 , height = 4 )

## 提取TP53，看Trunk的二次打击和Pivate的二次打击的差异
result$bialleic <- ifelse( result$CNV_type %in% c("LOH" , "Loss") , "Bialleic\n(Mutation with LOH or Loss)" , "Other" )
result2 <- result %>%
group_by(Gene_Symbol , TrunkUse , Type , bialleic , count_all ) %>%
summarize( count = sum(count) )
result2 <- data.frame(result2)

tmp <- c()
uniq_type <- c("All","IM + DGC","IM + IGC" )
for(typeN in uniq_type){
    print(typeN)

    tmp_res <- subset(result2 , Gene_Symbol == "TP53" & Type == typeN)
    p <- fisher.test(matrix(c(tmp_res$count) , ncol = 2 ))$p.value
    if( p < 0.001 ){
        p_text <- trans(p)
    }else{
        p_text <- paste0( "P == " , round(as.numeric(p) , 2) ) 
    }

    tmp_res$p <- ""
    tmp_res$p[1] <- p
    tmp_res$p_text <- ""
    tmp_res$p_text[1] <- p_text
    tmp <- rbind(tmp , data.frame(tmp_res))
}

col <- c(
    rgb(190,97,99,alpha=255,maxColorValue=255),
    "grey"
  )

names(col) <- c("Bialleic\n(Mutation with LOH or Loss)" , "Other")
tmp$ratio <- tmp$count/tmp$count_all
tmp$ratio[is.na(tmp$ratio)] <- 0
tmp$value_text <- round(tmp$ratio , 2) * 100
tmp$value_text <- ifelse( tmp$ratio == 0 , "" , tmp$value_text )

plot <- ggplot( data = tmp , aes( x = TrunkUse , y = ratio , fill = bialleic ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion (%)")+
    facet_grid(.~Type,space='free_x',scales='free_x') +
    theme(panel.grid = element_blank())+
    scale_fill_manual(values=col) +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 8,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold') ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Bialleic.TP53.trunk_private.pdf")
ggsave( out_name , plot , width = 7 , height = 4 )
